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Abstract. We address the issue of cosmological backreaction from non-linear structure for- 
mation by constructing an approximation for the time evolved metric of a dust dominated 
universe based on a gradient expansion. Our metric begins as a perturbation of a flat 
Friedmann-Robertson- Walker state described by a nearly scale invariant, Gaussian, power- 
law distribution, and evolves in time until non-linear structures have formed. After describing 
and attempting to control for certain complications in the implementation of this approach, 
this metric then forms a working model of the universe. We numerically calculate the evo- 
lution of the average scale factor in this model and hence the backreaction. We argue that, 
despite its limitations, this model is more realistic than previous models that have confronted 
the issue of backreaction. We find that the instantaneous effects of backreaction in this model 
could be as large as ~ 10% of the background. This suggests that a proper understanding 
of the cumulative effects of backreaction could be crucial for precision cosmology and any 
future exploration of the dark sector. 
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1 Introduction 

The issue of cosmological backreaction [1] has fueled a lively debate in the literature. The 
question is: can inhomogeneities in our universe backreact and affect the average dynamics 
within our causal horizon such that the observed acceleration can be attributed to their influ- 
ence? If the answer is positive then the mystery of the cosmological constant will have found 
a solution requiring no new physics but which will nevertheless demonstrate the subtlety 
of gravitational phenomena. Even if the answer is negative, backreaction will be operative 
at some level due to the non-linear nature of gravity and its effects may be visible in the 
new generation of cosmological observations. Either way, calculating its magnitude is an 
important cosmological question. 

The problem with assessing the magnitude of this backreaction lies with the complexities 
of the non-linear structures forming under gravity in the late universe. For example, using 
second order perturbation theory falls short of providing an answer since the effect is expected 
to emerge in the non-linear regime where perturbation theory breaks down. Attempts to 
model the non-linear structures involving toy models of voids and overdensities may be useful 
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for developing intuition but are nevertheless simplistic. N-body simulations imply that the 
effect is small but they are Newtonian with zero global backreaction by construction; the 
backreaction is a purely relativistic effect [2]. A strong argument for the smallness of the 
effect uses the fact that even with non-linear overdensities, the local metric perturbations 
and peculiar velocities are still much smaller than unity (away from black holes) and that a 
perturbative FRW framework should therefore hold. However, counterarguments have been 
put forward involving subtleties in the choice of background employed in these treatments. 
For references finding a a small backreaction see e.g. [3-8] while arguments for significant 
backreaction can be found in e.g. [9-14]. For recent articles and reviews on the subject with 
more extensive references see [15-26]. 

In this paper we discuss cosmological backreaction in a novel manner, with a fully 
relativistic framework for a universe that begins as a perturbed Priedmann- Robertson- Walker, 
= 1, CDM universe. To be precise, we employ a gradient expansion to express the metric 
as a series of terms with an increasing number of spatial gradients and coefficients which 
are functions of proper time. The initial conditions are the standard adiabatic and Gaussian 
post-inflationary primordial perturbations. We use the synchronous gauge: our coordinate 
lines comove with CDM particles and our time hypersurfaces are labeled by their proper time. 
Of course the gradient series has to be truncated and thus does not capture the developing 
non-linearities entirely realistically. However, we argue that even if truncated this series can 
still provide a well motivated model for the true geometry which can be made increasingly 
accurate in principle. Furthermore, it extends into the non-linear regime, describing the 
collapse of initial over-densities and the rarefaction of initial under-densities which go on to 
form the voids dominating the cosmological volume.^ 

Here is the outline of the rest of the paper: In the following section we obtain the series 
solution for the metric in a gradient expansion using a Hamilton-Jacobi formulation. This 
approach, first developed in e.g. [28] and [29], simplifies the calculation significantly compared 
to a more straightforward expansion of the standard Einstein equations. Our treatment 
follows a slightly different logical development. Then, in section 3, we apply these results 
to the backreaction problem by numerically evaluating the evolution of the average scale 
factor and the backreaction parameter Q. We find that, up to certain qualifications which 
we explain, backreaction leads to non-negligible deviations from the unperturbed background 
model. Our results indicate that the effect might be relevant, indeed crucial, for precision 
cosmology and the exploration of the dark sector. We summarize and discuss our findings in 
section 4 where future directions are also laid out. 

2 The inhomogeneous CDM Universe as a series in spatial gradients 

The study of the backreaction of cosmological inhomogeneities indicates that relevant effects, 
if any, will necessarily be in the relativistic and non-linear regime, in the realm beyond cos- 
mological perturbation theory. In this regime deviations from a homogeneous FRW Universe 
can be studied using an expansion in spatial gradients. Furthermore, calculations beyond 
the lowest quasi-homogeneous order can be drastically simplified via the application of such 

'^The synchronous gauge and a gradient expansion to study backreaction was first used in [27] to argue 
that a non-perturbative approach is really needed to settle the issue. In the present work we use a truncated 
gradient series but also consider it in the non-perturbative regime, beyond its apparent radius of convergence. 
As explained in the text we argue that this can be used as a model for the metric of the universe which 
captures the evolution of over- and under-dense regions. 
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a gradient expansion in the Hamilton- Jacobi formulation of gravitational dynamics. We 
present this formulation below and use it to obtain the first terms in the gradient series for 
the metric. We then compare with an exact spherically symmetric solution to gauge the 
accuracy of the approximation and develop some intuition. The section ends with a descrip- 
tion of how to apply this gradient expansion for evaluating the backreaction of cosmological 
inhomogeneities by controlling a number of complications that arise in such an application. 

2.1 Hamilton-Jacobi for CDM 

The Hamilton-Jacobi approach requires a Hamiltonian formulation which in turn requires an 
action. The action for gravity and a non-relativistic matter fluid (dust) can be written as 



(2.1) 



where x is 9- potential for the 4- velocity of the fluid, = —g'^^d^x^ ^iid /3, the energy 
density, acts as a Lagrange multiplier whose variation ensures that U^U^ = — 1. Variation 
wrt X gives the continuity equation V^(/?(9^x) = 0, while variation of the dust part of the 
action wrt g^^ gives the usual energy momentum tensor T^y = pUfJJy. Gravity is described 
by the standard Einstein-Hilbert term. 

Let us now use the ADM decomposition for the metric and develop a Hamiltonian 
formalism. The metric is written as 

goo = -N^ + hijN'N^ , goi = HjN^ , gij = -f^j , (2.2) 

with the inverse 

By defining the canonical momenta 



6jij 2k N 
5S 



(2.4) 

"jT = p^^l + r^d^xdjX, (2.5) 
with 

Eij = 1% - V(,7V,) , E = h'^Eij , (2.6) 
we can bring the action to the canonical form 

S = Jd'x (^TT^^ + TT^^^ -NU- NiU"^ (2.7) 

= (vr^.TT^^^ - y) - + vr^^yi+y^a^ (2.8) 



where 



and 



Ui = -2Vfc7rf + n^dix (2.9) 
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As is well known [30], the action (2.7) does not define a conventional Hamiltonian system 
since solutions to the equations of motion must set the "Hamiltonian" to zero, V( = 0, a fact 
which reflects the time reparametrization invariance of General Relativity. One can however 
obtain a conventional Hamiltonian formulation by using as a time parameter one of the scalar 
fields of the system [30]; in this case the obvious choice is to use the hypersurfaces of x 
time hypersurfaces, setting ^ = 1, = 1 and diX = 0. Thus, the metric takes the form 

ds^ = -dt^ + J^j{t,x)dx'dx^ (2.10) 

and the spatial coordinate lines comove with the matter. We then impose the energy con- 
straint l/( = and from this equation determine the (non-zero) tt^ which now plays the role 
of the Hamiltonian density: —ir^ =71. In particular 

U = 0^n = ^^ij^''' [liklji - ^7u7m) - , (2.11) 

and the action becomes 

S = j dtd^x - ^ + 2NiykTT^'^ . (2.12) 

In this form it defines a constrained Hamiltonian system where the canonical momentum tt*-' 
is constrained to be covariantly conserved 

Vfc7r'^^ = 0. (2.13) 

Let us now apply the Hamilton- Jacobi approach to the Hamiltonian system (2.12). 
Writing 

TT ■' = 

we obtain the Hamilton-Jacobi equation 



(2.14) 



which is a single partial differential equation for 5 as a functional of 'jij and a function of t. 
Once 5[t,7jj] is determined the metric can be obtained from 

Let us now turn to the remaining constraint (2.13). It will be automatically satisfied if 
S = f d^x.^ F{t,"fij) where J-{t,^ij) is a scalar function of the metric making S invariant 
under 3-D diffeomorphisms. Indeed, the variation of such a functional wrt the metric will yield 
a covariantly conserved tensor which will thus satisfy (2.13). In 3 dimensions all information 
about the spacetime curvature is contained in the Ricci tensor and T can be written as 
[28, 29] 

T = -2H{t) + J{t)R + Li{t)R^ + L2{t)R'^Rij + . . . , (2.17) 

a series in powers of the 3-D Ricci curvature involving an increasing number in gradients of 
7jj. H{t) will turn out to be the Hubble rate. Note that under our assumptions this form 
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is essentially unique. The Hamilton- Jacobi equation (2.15) can now be solved separately for 
terms with different number of gradients. This gives for the time dependent coefficients in 
(2.17) 



dH 
~dt 



dL 



dt 
dM 
dt 



dt 2 









L2H + 2 = 



e.t.c. The solution for H then determines all other functions: 

2 

3t ' 



H 



^=10*' 
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(2.18) 
(2.19) 
(2.20) 
(2.21) 

(2.22) 
(2.23) 
(2.24) 



Given the above, the metric can obtained by solving equation (2.16), which, up to terms 
with 4 spatial derivatives reads 

97, 



dt 



'-^ -= 2H-fij + J{Rj,j-4R,j) 
+Li (37yi?2 _ 8^^.^. + .^.) 

+L2{^^ijR^^ Rki — SRikR^j + ^ijR^^\k 

+4i?,'^|,fc + 4i2/|,fc - +4i?,/>) + . . . , (2.25) 

where the ellipsis denotes terms with 6 ore more spatial derivatives. This equation can now 
be solved iteratively. Setting an initial metric kij at time t\ and writing 



„(o) 



.(2) , Ai) 



we obtain at 0th order 



at 1st order 



and at second onder 



4/3 



+ 



(0) 



k 



„(1) 



^r-\ t- 
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(2.26) 
(2.27) 
(2.28) 



7. 



(2) 



81 



350 \t 



8/3 



tl 



-AR^'^Rim + -R-k'^ + ^-^^ ) 



lORRij + 17R iRij — -Rij-k ' + 8 



(2.29) 

where we have ignored terms that are subdominant when t ^ t\? An important point to 
note is that all hatted quantities are evaluated from the initial metric kij: R = R[kij]. 



Such terms ensure that hm 7^, = kij 



- 5 - 



Before proceeding, it is worth estimating the range of vahdity of the above gradient 
series approximation. In solving (2.25) we assumed that terms with more gradients are less 
important and that an iterative solution can be constructed. Prom (2.27) and (2.28) we can 
then estimate that the expansion should be valid for t < icon where 

Un-0(l-3)^. (2.30) 

The exact coefficient is of course dependent on the form of the initial 3- metric kij. After this 
time the contributions of successive terms to the metric containing more gradients are no 
longer perturbatively ordered and the iterative solution to (2.25) breaks down. However, we 
will argue below that this is enough to determine the non-linear collapse of over-dense regions 
and goes some way in describing the expansion of under-dense regions into the non-linear 
regime. 

Let us now apply the above to our universe. Assuming the standard inflationary initial 
conditions, the initial seed metric takes the form 

/ 10 



kij =[j-J + Y^(x) j (2.31) 

where $(x) is the primordial Newtonian potential and we have scaled the seed metric such 
that in the absence of perturbations the scale factor would be unity today. Note that the 
lowest order in the above expansion is simply the separate universe approximation where 
each spatial point scales as a FRW, CDM, universe with scale factor a cx 

In a flat, homogeneous, universe, to (the age of the universe) can be determined from 
the expansion rate, Hq = a/ a, or density, = l/(67rGto), observed today. When there is 
backreaction, assigning to is not straightforward. Ideally we wish to set to such that when 
t = to we obtain the average expansion rate and average density observed today. However, 
until performing the full numerical calculation we do not know what the resulting average 
values of these quantities will be. The problem is exacerbated by the transfer function used 
to obtain $(x) because it will also depend on the value of to. Therefore, altering to does 
more than just rescale time. 

We choose to define to through the Hubble parameter of the background, homogeneous, 
universe when it is extrapolated to t = to. For Hq^^^ = lOO/ikms"^ Mpc~^ we assign 

to = (2.32) 

Throughout, we take h = 0.4. It is true that in the observed universe indications are that 
h ~ 0.7 [31]. However this result is derived assuming a ACDM model where backreaction is 
assumed to be negligible. Our background model is a flat CDM universe; therefore if we take 
/i = 0.7 the density of that background universe today would be much larger than what is 
observed. If the backreaction is a small effect then the current average density of the universe 
will also be too large. The choice h = 0.37 corresponds to a background model that has the 
same density today as the best fit ACDM model [31], albeit a different expansion rate. If 
the backreaction is small then the average density in the model today will be close to the 
observed average density. This choice is of course arbitrary and any value for h in the range 
(0.37 — 0.7) would be equally appropriate. We have tested values throughout this range. 
The qualitative nature of our results and thus the conclusions we draw from them do not 
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change - in particular, while backreaction can be non-negligible in our model, it can never 
simultaneously reconcile Hq^^ and with observations. In future work, with A included, 
or with backreaction large enough to mimic A, it would be interesting to scan through all 
the relevant cosmological parameters, including h, to find the parameter set that best fit the 
observations. 

With kij as the initial metric, we have from (2.27) - (2.29) 



„{0) , ^,(1) , ^,(2) 



t \ 



4/3 



^ N 2/3 / i \ 



(2.33) 



where ^ 
^ 28 



19^,u^'\j - 12^, + ^Sij (($'\/)2 - $,/™^''")] (2.34) 
Given these approximations, the local density of matter can be obtained from 



p{t, x) = , I,, , ^ (2.35) 



'Det 



Note that if this expression is expanded to linear order one recovers the result of linear 
perturbation theory for the density contrast in the synchronous gauge. Equations (2.33) and 
(2.35) are the main results of the gradient expansion we will he using in this paper. 

The timescale for which the above approximation for the metric is accurate is estimated 

to be 

^=.3.4x— (2.36) 

in accordance to (2.30). We see that for describing an inhomogeneous universe which re- 
sembles ours for the whole of to) the metric can include perturbations down to about 
k = 0.3/1 Mpc~^ Of course, even shorter scales can be accurately described but for shorter 
times. It is then plausible that the gradient expansion can accurately follow the whole his- 
tory of our model universe at least containing inhomogeneities that have recently entered the 
non-linear regime. 

Unfortunately, to properly study backreaction we need to include perturbations over 
a broader range of scales and follow them for the whole of to- It might seem that the 
gradient expansion is of limited use for this and that all hope of studying backreaction in 
this framework is lost. However, as we will argue below, the situation is not completely 
irrecoverable. With certain modifications the above metric can still qualitatively capture the 
behaviour of collapsing structures and can also qualitatively capture the behaviour of regions 
of the universe that empty out and begin to expand faster than the background for the whole 
of to- 

When regions of the universe collapse, we can see that the collapse time roughly equals 
icon and thus for these regions and one can reliably follow the collapse with the above formal- 
ism. The particles in these regions gather in a volume that is effectively zero compared to the 
volume of the non-collapsed regions of the universe and once this has happened a collapsed 
region decouples form the rest of the universe. At these points, the precise volume or the full 
behaviour of the metric is irrelevant for calculating backreaction. To model the backreaction 
in the universe accurately, the metric needs only to correctly label which regions will collapse 
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and capture the behaviour of the expanding regions. We will discuss ways to achieve this 
below. 

To go some way towards addressing this issue we will now apply the above approxima- 
tions to a small initial spherical perturbation and compare the result to the exact solution. 
This will give at least a qualitative understanding of the behavior of the truncated metric 
and allow us to develop some intuition before considering the problem of backreaction. 

2.2 The spherical case 

In order to better assess the accuracy of the series expansions (2.26) - (2.29) and obtain 
some intuition we now compare with a simple model: a spherical concentration of dust 
representing an over-(under-)density. For simplicity we will assume that the over-(under- 
)density is uniform and thus characterized by a positive (negative) constant spatial curvature. 
The exact spatial metric inside the sphere takes the FRW form 



dx^ = a(r)^ 



2.37) 



1 — Kr^ 

where the dynamics of the scale factor a{t) are given by the Friedmann-Lemaitre equation 



d \ ^ 8ttG fo K 



3 a? 

whose solutions can be written in parametric form 



(2.38) 



aiu) = ^ (coshn — 1) , \/\k\t(u) = (sinhu — u) , (k < 0) (2.39) 

6 \k\ 6 \k\ 

aiu) = — (1 — cosn), \/\k\t{u) = — (u — sinu), (k > 0) (2.40) 

On the other hand, applying formulae (2.26) - (2.29) we obtain as approximations for the 
scale factor 




iM = (-| l/iTTTrKKol-l +~.'<f4\-] +^^^ (2-41) 



350 8 ° V'^o J 



where 



The comparison of the approximation to the exact solution is shown in figure 1. It 
is immediately apparent that the first term in the series already goes beyond linear theory 
in capturing the eventual collapse of a positively curved region (an initial overdensity in a 
perturbed FRW model). Additional terms improve the convergence to the exact behavior 
and predict the collapse time more accurately. This is in agreement with the expectations of 
the previous section that the series converges to the exact solution for times comparable to 
the collapse time. 

For regions with negative curvature (modeling initial underdensities) the first term is 
pretty close to the exact solution, reproducing the asymptotic a oc r behavior of an eventually 
empty region which is simply Minkowski space written in expanding coordinates. Although 
the additional terms improve the convergence initially, it is clear that at late times these terms 
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Figure 1. The radius for tlic ?' = shell of a spherically symmetric ball of dust with uniform density 
of positive (left) and negative (right) 3-curvature as a function of -\/Ht. We have set SttG = 3 and 
/o/k = 1. Larger values of /o/k, representing smaller curvatures, simply scale the above graphs to 
longer times and larger scale factors. The solid line is the exact solution, the dotted line is the 0th 
order approximation for the metric (flat FRW or "separate universe approximation"), the dashed 
line is the 1st order metric (with 2 spatial derivatives) and the dash-dotted line is the 2nd order 
metric (with 4 spatial derivatives). For positive curvature the addition of successive terms improves 
the convergence. For negative curvature, successive terms eventually deviate from the exact solution 
as the radius of convergence for the series is exceeded. However, a good description is obtained for 
relatively large values of -^/IkIt. 

deviate from the asymptotic behavior. Close examination reveals that this occurs at around 
the collapse time TcoI ~ -^j2 of the corresponding configuration with positive curvature, as 
expected in the previous section. However, the underdense regions are still reasonably well 
described up to times much longer than the collapse times of regions with corresponding 
positive curvatures: from (2.41) we see that at r ~ IOtcoI the solution deviates from the 
exact one by a few percent. 

The success of the gradient expansion for this simple case gives us confidence that 
our metric, eq.(2.33) can be used to model, at least qualitatively, the salient features of 
inhomogeneities in our universe that are relevant for backreaction: the collapse of overdense 
regions and the expansion of underdense ones significantly beyond the linear regime and in a 
manner more realistic than the models used so far. 

2.3 Potential complications for application to backreaction 

The previous subsection has indicated that the metric obtained form the gradient expansion 
can be a good model for studying backreaction. However, there are potential pitfalls when 
applied to more realistic distributions of matter that must be tackled before applying it to 
this problem. Firstly, in the real universe the metric at each point will not be spherically 
symmetric. Most importantly, there will be off-diagonal elements in the metric that are not 
present in the simple case of spherical symmetry. As we will see in section 3.2.4 these terms 
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can have a significant effect on tlie behaviour of the truncated metric. In section 3.2.4 we 
attempt to address this comphcation. 

Also, in the real universe, neighbouring regions interact. An initially underdense ex- 
panding region will not continue to expand as if it were isolated. Eventually it will expand 
far enough to hit other structures, which will necessarily affect the expansion. One would 
expect that as long as the series is convergent and with a sufficient number of terms retained, 
the metric in the gradient expansion should capture this behavior. However, for practical 
reasons we need to truncate the metric at a low order in gradients and eventually exceed the 
times for which the series converges. The fact that a low order truncated metric correctly 
describes the behavior in an isolated universe, as shown above, is not proof that it will also 
work well when there are interactions between neighboring regions. As we will see in section 
3.2.2, the expected behaviour of expanding regions in the real universe, when the regions 
interact, is not perfectly captured by the truncated metric. In section 3.2.2 we also attempt 
to address this complication. 

Finally, the synchronous gauge allows coordinate singularities where the determinant of 
the metric becomes and the density apparently diverges. These occur when the trajectories 
of particles, and consequently the coordinate lines they define, cross and caustics form. It is 
difficult to distinguish such a coordinate singularity from a region that is actually collapsing. 
One possible method is to calculate the Ricci scalar of the full 4-metric. At a coordinate 
singularity, the 4-Ricci scalar will not grow large, whereas in a true collapsing region it 
should at least grow larger than the background. However, even if regions could be correctly 
identified as coordinate singularities it would still be difficult to deal with them. Arguably, one 
could transfer into a more stable coordinate system, use that to travel through the singularity 
and then return to a new synchronous coordinate system. Remember, however that at these 
coordinate singularities particles are crossing each other's paths. There will inevitably be 
significant non- gravitational effects at these points. Even for weakly interacting dark matter 
particles, the single velocity fluid description will break down at these points. These effects 
are difficult to model, but they will certainly inhibit the expansion of the volume occupied 
by the crossing particles. 

To properly compare this model to the observed universe, this issue will need to be 
addressed more precisely.^ However, for our purposes, because we expect the dynamics of 
collapsed regions and regions with shell-crossings to be similar,^ we have chosen to label every 
region where the density diverges as a collapsed region. We expect that since any backreaction 
will be dominated by the expanding regions the accurate modeling of the contracting ones 
will not be critical. We explain in section 3.2.2 how we deal with the evolution of contracting 
regions in our model. 

2.4 Backreaction in the synchronous gauge 

Let us now briefly recall the definitions for backreaction and apply them to the synchronous 
gauge. Consider a fluid filled universe with the four-velocity of the fluid. The local 
expansion of the fluid lines is given by 

= u''.^ (2.43) 

•^A useful approach for future improvements might be an adaptation in this setting of the "adhesion 
approximation" described in [32]. 

*Even if their expansion is inhibited for different reasons. 
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while the stress tensor reads 



(^tiiy = ^ (uu;,! + u^-y) - {g^j,u + u^Uy) 9. (2.44) 
When one averages aver a domain D the average scale factor 

aDitf= I (fx^Z-yit,^) (2.45) 

JD 

obeys the Buchert acceleration equation [1] 

3^ = -47rG{p)D + Qd (2.46) 
an 

where 

Qd = 1 {{e')D - {0)1) - {a^ya^D (2.47) 
is the backreaction term. The brackets denote a spatial average 

{S)d = ,3^^ . (2.48) 

The above equation differs from the standard FRW one by the Q term, making it obvious 
that inhomogeneities affect the evolution of the average scale factor. In the synchronous 
gauge it reads 

Qd = \{{f'i^,)')D - - \{l''i^Jl^'ilk)D (2.49) 

As is implied by the notation above, the average quantities depend on the averaging domain 
and therefore the scale on which the averaging is performed. From now on we will mostly 
drop the the explicit reference to D but the inherent scale dependence for the quantities we 
refer to should be kept in mind. 

We can now use the expressions for the metric described above to estimate the amount 
of backreaction in our universe. We will directly evaluate the behavior of the average scale 
factor and cross-check with an explicit determination oi Qd. The evaluation requires the 
determinant of the metric and its inverse 7*-'. Even though we have an explicit form for 
7ij in eq.(2.33), writing an explicit expression for the determinant and the inverse is not 
convenient. Note that we cannot expand in a series since we are interested in the regime 
where the corrections from the gradient terms become order 1. It is however straightforward 
to approach the problem numerically as we describe in the following section. 



3 Numerical implementation 

To calculate the backreaction numerically, we set up an N"^ grid. We assign the initial 
gravitational potential $(x) on this grid to be a gaussian random field consistent with the 
known cosmology. Then, the evaluation of two spatial derivatives provides the matrices ^^ij 
and Bij in (2.33). Once these are obtained, expression (2.33) defines at each point on the grid 
a 3 X 3 symmetric matrix which provides an explicit time dependent model for the 3-metric. 
It is then straightforward to find the determinant 7 at each point and hence the average scale 
factor ao{t) by integrating over the points of the grid. Furthermore, the inverse 7*-' can also 
be calculated at each point and hence Q can be obtained. 
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3.1 Setting up the initial gravitational potential <I'(x) 

The initial metric perturbation is (approximately and sufficiently for this work) a gaussian 
random field. Since we'll be working in a box of size L we use a fourier series 

$(x)=^$ke*'^'' (3.1) 

k 

with periodic boundary conditions such that k = ^ {l,m,n), with l,m,n = 0, ±1, ±2, . . .. 
Reality of the field implies = $_k- Writing <l>k = -^k + ^-^ki the real and imaginary parts 
i?k and /k are gaussianly distributed random variables. The correlation function in a specific 
realization of the random field is given by 

^(r) = (<l>(x)<l>(x + r))v (3.2) 
where (. . .)v is a volume average. We then have 

e(r) = Y.{Rl + II) e^*^-^ (3.3) 

k 

The variance of and Ik will be {R^)e = (-^k)^ ~ ^"^fc' depending only on the magnitude 
of the wavevector. The ensemble average of the correlation function then reads 

(e(r))s = J^^rie^l^- (3.4) 
k 

Given a large number of modes and the ergodic assumption, the correlation function measured 
in our universe £,M{f) can be interpreted as (■^(r)) e and thus cr^ is determined by the measured 
power-spectrum. If our box is large enough such that ^ — t- j^tj J (i^k we have 

4 = 1^^^ 1—] T{kf (3.5) 



25 \^ Y ^piv / 

where T{k) is the transfer function that evolves the primordial gravitational potential through 
the radiation dominated era. Throughout, we use the transfer function of [33] with the 
baryonic correction presented in [34]. For the primordial spectrum we take, at the pivot scale 
^piv = 0.002 Mpc, the spectral index to be ns = 0.963 and the amplitude = 2.42 x 10"^ 
as measured by WMAP [31]. As inputs to the transfer function we also take h = 0.4 and 
= 0.17.^ We choose this value of Qh because it preserves the ratio between ^If, and ilm 
observed in the real universe. 

To create a given realization one has to draw i?k and Ik of each fourier mode in the box 
from a random process with a gaussian probability distribution 



/2 

P(Xk)dXk = exp 
0"fcV27r 



(3.6) 



The probability dP to obtain a given value of the complex <I>k is 

dP = P{Ri,)P{I^,)dR^dI^. (3.7) 



'Each of these are quantities of the background metric at t = to, not the full backreacting CDM universe. 
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Changing variables to <I>k = e*^'' we get 



dP = — ^ exp 



dA^ ^ . (3.8) 



We see that the amphtude is described by a Rayleigh distribution and the phase is uniformly 
distributed in the range [0, 2ii). Using this, we can also set up a realization of the gravitational 
potential and its derivatives from 

$(x) = 2 ^ Ak cos (k • X + 6'k) , (3.9) 

kGuhs 

and 

$ •j.(x) = -2 ^ hkjAi, cos (k • X + (9k) , (3.10) 

kfEuhs 

This provides all the necessary information for the evaluation of the metric at each point for 
all times using eq.(2.33). 

Note that, for the derivative $^ij(x), the sum is only very weakly convergent. It is 
instructive to count, in the continuum and large k limit, the factors of k involved in the 
variance of ^^^^(x). From a'j, there are factors of k~^ , T{k)'^ and k^"~^, from the second 
order derivative there are two factors of /c^, and from the integral measure there is a factor 
of k^. In the large k limit, T[k) — )■ ln(A;)/A;^. Therefore, the value of the partial sum, up 
to a value of k, scales as ln(A;)^ x While in principle the red-tilted power law k^"~^ 

will eventually fall more quickly than the growing ln(A;)^ term, this will only happen at very 
small scales (because In^ — 1| <C 1). This becomes important in the following section. 

3.2 Dealing with the complications and presenting the backreaction 

Ideally we would like to use eq.(3.10) to generate <I>^ij(x), then substitute this into eq.(2.33) 
and finally use eqs.(2.45) and (2.49) to evaluate the average scale factor ud and the magnitude 
of the backreaction, Q, in our choice of CDM universe. However there are a number of 
complications that arise before this can be straightforwardly achieved. Some of these were 
flagposted in section 2.3 and relate to the truncation of the metric or the extrapolation to 
long times, the use of the synchronous gauge and the almost divergent nature of the sum in 
eq.(3.10). 

The gradient formalism presented so far involves no assumptions about the universe 
except that it is dominated by dust, was once flat and homogeneous and is well described by 
general relativity. Every new feature that we must now add to deal with the complications 
arising adds extra assumptions about the universe. Therefore we are not able to calculate 
the backreaction from all scales in a dust filled universe with certainty. However, we can 
provide a plausible candidate model to describe the backreaction. The backreaction in the 
CDM version of this model is not large enough to mimic A but would be large enough to be 
measurable. If this model were extended to include A it could be tested against observations. 
It would be interesting to see if it provides a better fit than standard ACDM for a wider set 
of data. 

We discuss in section 4.1 various ways in which this model can be improved and refined 
beyond just the obvious inclusion of A. Nevertheless, we argue that even in its current form, 
the model resembles the actual universe more closely than earlier models in the literature 
that deal with backreaction. Once the initial conditions are set up, each grid point evolves 
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independently from the others. However, we show in section 3.2.6 that there is a depth of 
structure in the evolving box of grid points. This is because the initial conditions of the 
metric at each gridpoint are not independent. 

In the rest of this section we proceed to discuss the complications encountered when 
applying the gradient formalism to the problem, our methods for dealing with them and the 
corresponding results for the backreaction. 

3.2.1 The almost UV divergence 

The sum in eq.(3.10) converges as ~ ln(A;)/A;(^~"=)/^. However ln(fc)/fc(-^~""'^/^ does not begin 
to decrease until k ^ 1. It is highly likely that other natural cutoffs such as the free- 
streaming length of dark matter will regulate this sum long before it converges. Therefore 
it is not tractable (or even physically motivated) to calculate the value of this sum in the 
limit A; — )• oo. More importantly however, structure growth at any given scale will freeze 
out after some time. No kpc sized scales are still evolving gravitationally today. They have 
either collapsed into bound structures or are part of much larger (in synchronous coordinates) 
expanding regions. 

There are also practical limitations. Firstly, a box of A^^ grid points with a grid spacing 
that corresponds to the natural cutoff of eq.(3.10) will either need to be very small, or have 
very many points in it. We obviously cannot have an arbitrarily large number of grid-points 
without an arbitrarily powerful computer. If instead we choose to examine a small box then 
we lose information about correlations over distances greater than the box size. The model 
therefore loses the ability to calculate the backreaction that arises from structure growth over 
those scales. Secondly, we saw above that under-densities on short scales cannot be accurately 
followed by the gradient expansion for a arbitrarily long times unless extra assumptions are 
introduced made about their late time behaviour. Thus, even if we could include all scales 
from the very long to the very short in eq.(3.10) it would be unreasonable to expect that our 
approximations would make sense. 

We are interested in the scales that are evolving most significantly today. Therefore, 
we choose to examine boxes where the grid-spacing corresponds to these scales. The near 
divergent nature of eq.(3.10) means that it will then be these scales that dominate the value 
of $^ij(x) and thus the behaviour of the evolving metric. By varying the number of grid 
points we then vary the size of the box. For a box with grid spacing R and number of points 
Nx its size is L = RN^- We then sum eq.(3.10) from the value /cmin = (27r)/L up to the value 

fcmax = {27t)/R. 

Although this is necessary to have a working model, it does have the drawback that all 
evolution of the universe at smaller scales is removed from the calculation. Our model sees 
the universe as a collection of synchronous "blobs". These blobs expand and contract, but 
any internal structure of each blob is ignored. Similarly, as argued above, any backreaction at 
much larger scales than the grid spacing is not seen because of the finite box size. We cannot 
therefore capture all the backreaction in the universe over all times with our model. What we 
hopefully do capture is the backreaction that occurs due to the scales we are examining. In 
our plots we choose to focus on approximately the scales that have the largest backreaction 
from t ~ to/2 until t = tQ. This corresponds to ii = lMpc//i. In figure 5 we show the 
backreaction over a range of grid spacings. In section 4.1 we discuss possible methods to 
construct a model in the future that accounts for backreaction over all scales. 
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3.2.2 Fixing the metric 

Over-densities: If we let the metric evolve without any alterations then in all of the 
collapsing regions (i.e. grid points with increasing density) the determinant of the metric 
will eventually change sign. When this happens, both the local density and local volume will 
become imaginary. Clearly this is both unwanted and unphysical. In the real universe when 
a region collapses it virialises and its growth is subsequently frozen. We choose to model 
this by freezing the metric of any grid point when the ratio between the density at this grid 
point and the background density has crossed a certain threshold, prat-^ This corresponds to 
a condition on the determinant of the metric, ■~f{t,x) at each grid point of 



For every grid-point, x*, where this condition is satisfied at some time t*{x*), we define 
7jj(t,x*) = 7jj(t*,x*) for all t > t*. By extension, for the purposes of calculating Q, we also 
define 7jj(t,x*) = for all t > t*. We find that once p^nt is set to a sufficiently large value 
the precise value of pra.t has little effect on the results for Q and an- This is not surprising 
given that, independent of the value of the threshold, the asymptotic state for any of these 
grid-points is a point of zero volume relative to the volume of all of the expanding grid-points. 

Under-densities: We showed in section 2.2 that the truncated metric describes the evo- 
lution of an isolated, under-dense, spherical region well. Unfortunately, the universe is not 
made up of a collection of truly isolated points. If the truncated metric is left unaltered 
then expanding regions will continue to expand faster than the background, indefinitely. 
When keeping just one gradient term in the truncated metric, the asymptotic state of each 
under-dense point is equivalent to an empty, open universe. When keeping two terms the 
asymptotic state expands even faster. This behaviour results from the series not converging 
anymore on such timescales. In the real universe the expanding regions will eventually hit 
walls and filaments that will impede (or perhaps even promote) their expansion. At this 
point the region corresponding to a given grid-point will continue evolving as a part of a 
bigger over-all region. That is, larger scales will have become relevant for studying structure 
growth. 

We do not know what will happen to each grid-point when these larger scales become 
relevant. The only way to know is to "zoom out" and recalculate the universe in a larger 
box, with a larger grid-spacing. Unfortunately, different under-dense grid-points will expand 
at different rates. Therefore, they will expand into neighbouring regions at different times. 
So, we cannot just evolve the metric until larger scales become relevant and then stop. Some 
regions will reach this point much earlier than others and the most extreme regions will reach 
this point very early on. At the very least, we need a way to deal with the most extreme 
regions, but ideally we also want to know at all times what magnitude of backreaction has 
been generated by the scale being examined. This includes times long after other scales have 
also become relevant. 

We attempt to deal with the expanding regions by taking any gridpoint that has over- 
expanded and fixing it to evolve with the background expansion. We tag a grid-point as 

^See section 2.3 for a discussion concerning the potential misidentification of coordinate singularities as 
collapsed regions. 




(3.11) 
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over-expanded when the ratio of its volume to the comoving background volume at that 
gridpoint has exceeded a certain threshold, ^rat- That is, when 

7(t,x)><7?at(^) • (3.12) 

For every gridpoint, x*, that satisfies this condition at some time, t*(x*), we fix the metric 
to have the following behaviour 



7ij = , i / j 

. _ 4 

la — 3^^** 

iij = 0,i^j (3.13) 

for all t ^ t* . Note, 7 = 7(t*,x*). If we impose this condition immediately at t = t* the 
discontinuous change in 7 causes numerical artifacts to appear in Q. Therefore we actually 
impose this transition smoothly in a manner described in appendix B, with eq.(3.13) only 
the asymptotic final state of any over-expanded gridpoint. This transition can be made 
arbitrarily fast; however see appendix B for a discussion. Unlike the situation for collapsed 
regions, our results are highly sensitive to the value of g'rat- This is also not surprising given 
that it is the expanding regions that come to dominate the volume of the universe. 

The motivation for this particular method of fixing expanding regions is that once a 
region has crossed this threshold, any backreaction that occurs is due to evolution on larger 
scales. Therefore the backreaction caused at this gridpoint by the scale of interest is, arguably, 
zero. By fixing the grid point to evolve with the background we force there to be no additional 
backreaction from that region. 

Note that even if the complete solution for the metric could be obtained it would still 
not be an entirely realistic description of the non-linear regime. Indeed the contraction of 
overdensities is eventually halted by further baryonic physics and structures virialize (un- 
less they are dense enough to form black holes) while underdensities are rarefied up to the 
point when the shells surrounding them cannot expand further due to the interaction with 
surrounding overdensities. Although apparently artificial, this treatment of non- linear scales 
should in principle capture the effect of the actual baryonic physics preventing gravitational 
collapse and leading to virialized structures which decouple from the background evolution. 

3.2.3 Two definitions of Q and a 

We want to calculate the values oi and Q in our model universe. The definitions for Q 
provided in eqs.(2.46) and (2.49) assume a universe composed only of dust. Ostensibly, this 
is true for the model. However, as soon as we fix expanding regions to the background by the 
"hand of God" there is a risk that the metric will behave in a way that no longer describes 
a dust only universe. If the universe behaves in a non-dusty manner then we expect that 
the Q calculated from the two equations will differ. This is because there will be an effective 
pressure term that should be present in eq.(2.46). 

We do not want our somewhat ad hoc fixing of the expanding regions to be interpreted 
as due to the effects of backreaction. Therefore, we calculate Q in two manners.^ The first 



''For a related Newtonian analysis see [35] 
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method, which we still label as Q, we calculate directly from eq.(2.49). The second method 
involves numerically differentiating aoit) twice to obtain aoit) and using this to define an 
effective Q{aD) from eq.(2.46). Correspondingly, we also calculate ao in two ways. The first 
is directly from the definition, eq.(2.45). The second method involves numerically integrating 
eq.(2.46) using Q as calculated from eq.(2.49). 

It is the calculation of Q{aD) that is most responsible for our need to impose the fix of 
section 3.2.2 in the smooth manner described in appendix B. To calculate Q{a£,) we require 
qd- If the evolution of 7 at any point is not smooth, then the evolution of a/) is not either. 
Thus, d£) is not continuous and 'do is not well-defined. Even with the smooth fix applied to 
the expanding regions Qiao) exhibits a large degree of numerical noise. In figures 3 and 4 we 
first smooth ao over a small range of time before taking the numerical derivative to obtain 
a£) and do- 
lt could be argued that the best method to fix jij in the expanding regions is the one 
for which both methods of Q and both methods of a/j align most closely. 

3.2.4 Backreaction with artificial local isotropy 

The determinant of the metric tensor is given by 

7 = 711722733 + 2712723713 - 712733 - 723711 - 713722 • (3.14) 
To first order in the expansion, the metric tensor itself is given by (see eq.(2.33)) 




(3.15) 



is distributed symmetrically around zero. Therefore, if we ignore correlations, the first 
term in 7 is symmetrically distributed around (t/to)^- The second term is symmetrically 
distributed around zero. However, is always positive. Therefore, the third, fourth and 
fifth terms in 7 will not be symmetrically distributed around zero. The net result is that 7 
is not symmetrically distributed around (t/to)^, but is strongly skewed towards 7 < (t/to)^. 

If we naively apply the formalism as it has been developed up to this point, then our 
model universe will consist almost entirely of collapsed regions. It is tempting to interpret 
this result as physically motivated; however it is not. It is an artifact of the truncation of the 
metric. While we have been consistent in keeping the same number of gradient terms at each 
order of the expansion in the metric, this is not true for the metric's determinant. Even when 
we include just one gradient term in the metric, the metric's determinant has terms up to 
cubic order in the gradient. It might then be tempting to also truncate the determinant at a 
certain order of the gradients; however this makes the model itself inconsistent. We may have 
truncated the metric, but the terms in the metric are not perturbatively ordered; therefore 
removing terms from the metric's determinant would mean that our model no longer satisfies 
general relativity. That is, we would not be using the correct volume element for the model's 
metric tensor when calculating Q and ao- 

We solve this, the last of the model's complications, in two ways.^ In this section, we 
consider a simplification of the metric such that = <I>^ii and = for all i 7^ j.^ 

^The second method is introduced in section 3.2.5. 

^Note that we do not leave "I>,22 and $.33 as their true values because the non-zero correlation between 
and also skews 7 towards collapsed regions. 
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This amounts to considering every grid-point in our model as a spherically symmetric region. 
While this choice of initial condition is less well motivated than using the full gradient 
once the initial condition is set up it behaves entirely consistently. Moreover, it provides a 
good point of comparison. We know from section 2.2 that the truncated metric describes an 
isolated spherical region well. Therefore, up until the point when larger scales should become 
relevant, we can be confident that the metric will describe the evolution of these initial 
conditions accurately. Although much of the information about the correlations between 
different grid-points is removed because of this simplification, we show in section 3.2.6 that 
not all of it is. 

Note that even in this simplified model, the distribution of itself matches that 
expected from a near scale-invariant, Gaussian, power-law primordial spectrum. Therefore, 
while information about correlations between points is lost, the points themselves should 
have the correct distribution. 



Figure 2: We have plotted 



to J 



2/3 



^« = ./.2/3 (3-16) 

to j 

for this model in figure 2. The grid spacing used for these curves is ii = lMpc//i and the 
number of grid points is A'^ = 30^.^'' We have used both the definitions of presented in 
section 3.2.3 and have shown results when the metric is truncated to include one gradient 
term or two. Finally, we have shown the results when we vary the threshold, (^rati described 
in section 3.2.2. As is expected from the earlier results in section 2.2, the order of truncation 
does not have a significant effect on the results. For the case where there is no fix at all on the 
expanding regions (upper panel) all four methods to calculate ao give similar results. While 
this result is unphysical because the expanding regions will not really expand indefinitely 
without interaction, this result reassures us that at the very least we have the mathematics 
of our model under control. 

The curves that result when we do fix the expanding regions to the background are 
also what we expect. After sufficient time has passed, the curves that use the actual average 
scale factor, a^, stop increasing. This time corresponds to the point when the bulk of points 
in the grid have either collapsed or been fixed. At this point any backreaction at this scale 
ends and Aa becomes constant. When the threshold is increased, this turnover occurs later. 
Also, due the non-dusty method of applying this fix, the curves that are obtained from Q 
do not see the effective pressure and continue increasing. Note that even once Q — )• 0, a{Q) 
continues slowly increasing. This is because a continuous pressure is required to keep the 
expanding regions from expanding faster than the original background model. 

We do not know what the correct g^at value is for our universe. In an extension of our 
model that includes A, the ao curves in figure 2 would serve as a prediction for the magnitude 
of backreaction expected in the universe as a function of g^at ■ 

Figure 3: We have plotted 

= (3-17) 



^"This grid spacing was used because, as can be seen in figure 5, it is approximately this grid spacing that 
gives the maximum rate of backreaction in the late universe. 
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Figure 2. The relative difference between the average scale factor, au, and the scale factor of the 
background model, [t/tof'/^^ is plotted against time for the model presented in section 3.2.4. This 
model consists of a grid of spherically symmetric regions, where = •I'.ii- $,ii is generated from a 
nearly scale-invariant, Gaussian, power-law, primordial spectrum. All three plots use the unsquared, 
metric in eq.(2.33). The metric has either been truncated after including one gradient term (i.e. 
or two. The upper panel has no fix on expanding regions. The lower two panels do, with thresholds 
of giat = 2 and 4 respectively (see section 3.2.2 for details). The ao curves use the actual average 
scale factor defined in eq.(2.45). The a{Q) curves use an average scale factor reconstructed from Q 
using eq.(2.46) (see section 3.2.3 for details). The grid spacing used in this plot is i? = 1 Mpc//i and 
the number of grid points is A^^ = 30^ (see section 3.2.1). to is the time today. 



for Q and Qiao) in figure 3 (see [36] for a definition of ilg). The model is identical to that 
used in figure 2. The left panel has no fix on the expanding regions and the right panel uses 
9rat = 4. We see the same behaviour in this figure as we did in figure 2. When there is 
no fix, all four methods for calculating Q roughly coincide. When we do fix the expanding 
regions to the background, Q is roughly unaffected. However Q{a£)) changes significantly 
due to the effective pressure this fixing produces. The magnitude of "pressure" being applied 
could in fact be calculated from the asymptotic value of ^Qi^ao]- Note that in all examples the 
maximum value of Q.q is > 0.1. Therefore, at least in this model, it appears that backreaction 
is large enough to have a measurable impact on dark energy constraints. Note here that even 
after od has formed a plateau, Q{a£)) is negative, this is because of the pressure required to 
keep the expanding regions expanding at the same rate as the background. 
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Figure 3. Q,q is plotted against time for the model presented in section 3.2.4. The grid used is 
identical to the one used in figure 2. The left panel has no fix on expanding regions. The right panel 
does, with a threshold of ^rat = 4. The Q curves use the actual value for Q defined in eq.(2.49). The 
Q{a£i) curves use a Q reconstructed from ao using eq.(2.46) (see section 3.2.3 for details), to is the 
time today. 



3.2.5 Backreaction with the squared metric 

Our initial motivation in using the gradient expansion was to construct a numerical model 
universe where the correlations between separate regions were taken into account. We de- 
viated quite significantly from this motivation in the previous section by removing a lot of 
information from the initial metric. This was necessary because if we used all of the com- 
ponents in ^^ij the determinant of the metric 7 was unphysically skewed towards collapsing 
regions. We explained in that section that this occurred because one of the terms in the ex- 
pression for 7 was unphysicallly skewed towards negative values as a result of the truncation 
of the metric. 

In [28] the authors performed a trick to avoid the determinant of this truncated metric 
crossing zero. We explain this trick and write the metric that it produces in appendix A. We 
will call the derived metric the "squared metric" because its derivation involves taking the 
square of the metric. The resulting metric contains at each order of the gradient expansion 
some extra higher order terms to ensure the positivity of the determinant. In [28] it was found 
that this squared metric actually improved the match between the gradient expansion and 
the exact solution they were studying (the Szekeres solution). Unfortunately, as we show 
in appendix A, the squared metric matches our test case, the expanding and contracting 
spherical regions, worse. 

Nevertheless, even at the first order of the expansion, the squared metric includes terms 
that involve the square of the gradient. Therefore the squared metric behaves much less 
pathologically. This is because, as well as the terms that are skewed towards negative values 
in eq.(3.14), there will now include many other terms that are skewed positive. Whether a 
region collapses or gets fixed as an expanding region will depends on which terms are larger 
in that region. The non-pathological behavour of the squared metric allows us to use all 
of the components of in the metric, thereby keeping all the information concerning the 
correlations between grid points. This also makes each individual grid point more realistic 
as it will not be isotropic and will expand and contract at different rates along different axes, 
which is exactly what happens in the real universe. 

Fixing the squared metric: There is a much more diverse range of behaviour with the 
squared metric. Some regions first expand, then contract, then expand again a number of 
times. For the one-gradient truncated metric, the asymptotic behaviour is however the same. 



- 20 - 



That is, some regions end up expanding indefinitely and others collapse to a point. We 
fix the expansion for this metric in exactly the same manner as before. However, for the 
two-gradient truncated metric, the asymptotic behaviour is different. No single region in the 
grid ever reaches a point of asymptotic expansion. This would be great if it were not for the 
fact that instead, all regions eventually collapse. This is clearly just as unphysical as every 
expanding region in the universe expanding indefinitely. 

We still wish to fix the expansion of any over-expanded regions to the background. 
However we choose to take advantage of the extra complexity afforded to us by this metric. 
Given that every initially expanding region will eventually collapse, we fix the expanding 
regions and freeze the collapsing regions in a two step process. Firstly, exactly as before, 
we set a density threshold for collapsed regions, /Orat; and a volume threshold for expanding 
regions, g^at- Once any gridpoint satisfies either the condition in eq.(3.11) or (3.12) we peg 
that region as either a collapsing region or an expanding region respectively. However, we 
do not yet alter its evolution in any way. Whenever the density at a collapsing grid point 
stops expanding more slowly than the background and begins to expand faster than the 
background, we freeze it. Equivalently, whenever an expanding grid point turns over and 
begins to expand more slowly than the background, we fix it to the background. That is for 
a region labelled as collapsing we freeze the metric when 

2f > (3-18) 

is first satisfied. And for a region labelled as expanding we fix the metric to the background 
when 

is first satisfied. Note that these expressions should be understood to be evaluated in the 
limit that At — t- 0. Note again that once these points are reached, the method we use to 
freeze and fix the metric is identical to that described in section 3.2.2 and appendix B. 

Figure 4: We have plotted against time for the model universe generated using the 
squared metric. For these curves we have used /?rat = 3 and g^^t = 3. R = lMpc//i as 
usual. What is most interesting about these curves is that, for the truncated metric with 
two gradient terms, Q{aD) and Q follow each other quite closely. This is a very nice result. 
The same is not also true for the truncated metric with one gradient term; however this is 
not surprising. Remember, the one gradient metric is fixed in exactly the same way as it 
was for the curves in figure 3. As a result the same type of effective pressure will be present 
causing these two curves to deviate. The fact that the two-gradient curves match indicates 
that in that model of the universe, the universe remains dust dominated, even when we fix 
the metric and there is substantial backreaction. 

The blue (solid) and red (dashed) lines therefore correspond to our best estimate of 
the backreaction in a true, statistically homogeneous, CDM universe that began fiat. In 
this model, the backreaction generated at a given scale is first positive and then negative. 
These periods correspond to an initial over-expansion of under-dense regions, followed by an 
inhibition of these regions' growth which occurs initially as a natural result of the metric 



^^Though note that if left alone, these "collapsed" regions would turn around and start to expand again 
which is an artifact of the squaring of the metric. This never happens in our model because we always freeze 
the evolution of any collapsed region. 
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Figure 4. ftg is plotted against time for the model presented in section 3.2.5. This model uses all the 
components of which is generated from a nearly scale-invariant, Gaussian, power-law, primordial 
spectrum. The metric used is the "squared metric" in eqs.(A.l) and (A. 2). The metric has either 
been truncated after including one gradient term (i.e. or two. The curves all use thresholds of 

<?iat = 3 and prat = 3 for the fixing and freezing of expanding and collapsing regions respectively (see 
sections 3.2.2 and 3.2.5 for details). The Q curves use the actual value for Q defined in eq.(2.49). 
The Q{aD) curves use a Q reconstructed from using eq.(2.46) (see section 3.2.3 for details). The 
grid spacing used in this plot is i? = 1 Mpc//i and the number of grid points is = 30^ (see section 
3.2.1). to is the time today. 



and subsequently due to our fixing the expansion to the background. The asymptotic value 
for both Q and Qiao) is approximately zero, indicating that eventually the backreaction 
generated at a given scale dies away. 

The magnitude of can be as large as ~ 0.2. Using different grid spacings gives 
qualitatively identical results, except for a rescaling of the time axis. Altering p^at and g^at 
can have a significant effect on l^igl; however it maximum value is always at least ~ 0.1. 

Figure 5: In a universe with significant backreaction the averaged Friedman equation also 
contains a term from the averaged spatial curvature (R) - see e.g. [36].^^ This term depends 
on the cumulative effects of backreaction; therefore it can be quite large and can contribute 
at least as much as Q to the measured effects of backreaction today. In this section, we 
display the total, cumulative backreaction. 

The total contribution of backreaction to the average expansion rate can be encapsulated 
in one term, defined through Qx = ^q + ^r- Note that Qr is defined in the same way as 
0,Q using the average expansion rate. Unfortunately, to calculate {R) at t > ti it is necessary 
to calculate second order spatial derivatives of the metric. This would involve calculating 
fourth order derivatives of the gravitational potential term that we use to generate the initial 
condition for the metric. We chose to circumvent this difficulty by directly comparing the 
density of matter to the expansion rate. This is possible because, in an initially fiat CDM 
universe, 

Qx{t) = 1 - nm{t) (3.20) 

^^Thanks to Thomas Buchert for pointing this out. 
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Figure 5. The cumulative backreaction J7x is plotted against redshift for the model presented in 
section 3.2.5. From right to left the solid blue curves correspond to grid spacings of i? = 0.01 Mpc/h, 
O.lMpc/h, IMpc/ft,, lOMpc/ft, and 100Mpc//i (in each case = 20). AU of the curves use the 
squared metric, truncated after including two gradient terms. The curves all use thresholds of (7,at = 3 
and Prat = 3 for the fixing and freezing of expanding and collapsing regions respectively (see sections 
3.2.2 and 3.2.5 for details). The dotted red line is the expectation in a ACDM model with no 
backreaction and fl\ = 0.7 today. 

where r2m(i) is defined by the average density and average expansion rate at time t. 

In figure 5 we plot this against "average" redshift, defined through = 1/(1 + z). 
We also show in this figure the efi'ects of changing the grid spacing. Each curve in this plot 
uses the same values of prat and gmt and is constructed using the squared metric truncated 
after including two gradient terms. It is clear that a smaller grid spacing induces backreaction 
at earlier times. One interesting feature of this plot is that while scales of approximately 
1 — 10 Mpc//i are showing relatively large backreaction today, when averaged over 100 Mpc//i 
scales the universe is clearly still very homogeneous with small backreaction. From the 
arguments made in section 2.1 we expect that up to t/tQ = 1 (z = 0) the region below the 
10 h-^ Mpc line is a realistic expectation for the magnitude of backreaction from such scales 
alone. However, this plot does indicate that a statistically homogenous universe does not 
preclude the possibility of significant backreaction arising from structure growth at smaller 
scales. For reference, the dotted red line is the expected curve in a ACDM model with no 
backreaction and = 0.7 today. 

The fact that a range of scales that vary by more than an order of magnitude can all 
be backreacting at the same time indicates that our use of one fixed grid to assess the total 
backreaction in our model universe is lacking. If a truncated gradient expansion is to be used, 
then to properly quantify the total backreaction occurring in the universe at any given time 
it will not be sufficient to focus only on a small range of scales. How to tackle this problem 
is an open question, although we suggest some possibilities in section 4.1. 
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Figure 6. This figure depicts the volume at each gridpoint in a two-dimensional slice of our iV^ grid 
at t = to- The upper left panel is an actual slice from our model, the other three panels are the result 
of redistributing each gridpoint to a new gridpoint using a completely random permutation. The axes 
show the synchronous comoving coordinates of each gridpoint, which therefore do not define physical 
distances at tg. Instead, the volume of each gridpoint is depicted by colour. The white regions are 
made of collapsed gridpoints with effectively zero volume. The blue and red (dark) regions are all 
expanding regions. The more red a gridpoint is the bigger the volume. This figure was produced using 
the model described in section 3.2.5. The list of parameters used for this figure appears in section 
3.2.6. 

3.2.6 What our model universes look like 

All of the figures we have presented so far have quantified in some way the magnitude of 
backreaction in our models; however they do not show very well what the model universe 
actually looks like. The initial conditions of the model include correlations between grid 
points, but the subsequent behaviour of each gridpoint is completely independent. The 
gradient expansion is an iterative solution to the full equations of general relativity, so a 
complete summation of this expansion should still describe the behaviour of each gridpoint 
exactly. However it is an important question to ask whether our truncated expansion retains 
information about the correlation between points, even after the terms in the expansion 
become large. 

In figures 6 and 7 we depict the volume of a two-dimensional slice of our grid for the 
models of sections 3.2.5 and 3.2.4 respectively. The gridpoints are defined by synchronous 
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comoving coordinates therefore the axes in these figures and the width of each region do not 
depict the area of each region. Instead, we depict the volume of each gridpoint using colour. 
The white regions consist of all the collapsed gridpoints. The blue and red regions are the 
gridpoints that have been labelled as expanding. The more red a coloured region is the bigger 
its volume is. For each figure we use a gridspacing of -R = 1 Mpc/h and = 50 gridpoints 
in each direction. We also use pmi = 5rat = 3 for each figure. The volume of each gridpoint 
is defined simply as the value of at that gridpoint. Finally, both figures are shown at 
t = to. 

In figure 6 we show a slice of the grid for the model in section 3.2.5. This is shown in the 
upper left panel of this figure. In the other three panels we have redistributed each grid point 
using a random permutation. It is reasonably clear from this figure that the structures in 
our model universe look different to structures in a completely random collection of over and 
under-dense regions. The collapsed regions clump together more, as do the over-expanded 
regions. It is even arguable that a few filaments of collapsed regions can be made out in the 
slice from the real model. 

In figure 7 we show a slice of the grid for the model in section 3.2.4. There is one 
striking feature of this figure that is immediately obvious, even to the human eye. There is 
a very clear preferred axis in this model. Almost all of the collapsed regions and expanding 
regions align along this axis. This is actually not surprising and in fact shows very clearly 
that correlations in our model do survive to late times. In constructing this model we set 
= ^,11- Therefore we lost all information about correlations in the X2 and X3 directions, 
but not in the xi direction. It is the xi axis that the structures line up along. One other 
feature of this figure is the greater uniformity in the volume of the expanding regions. This 
is also not surprising given how the expanding regions are fixed to the background in this 
model. In the model of section 3.2.5 we first label a gridpoint as expanding and then fix 
it only when its expansion slows to be the same as the background expansion. Different 
gridpoints will turn around when they have different volumes. In the model of section 3.2.4 
(used in this figure) the expanding regions are fixed as soon as they cross a certain threshold. 
The result of this and the method we use to do the fixing is that every region asymptotes 
to the same volume. The small differences that are seen in figure 7 result from some regions 
not having yet reached the asymptotic metric. 

The fact that non-random structures are clearly visible in our models suggests that the 
of use a truncated gradient expansion to describe the universe is not completely without 
merit. At the very least, the existence of these correlated structures is one step better 
than earlier models that are made up of completely independent regions. Our models lack 
the sophistication of the best Newtonian N-body simulations. However backreaction is a 
purely relativistic effect. Therefore our approach could at present be considered as the most 
sophisticated one. We do not require any special symmetries of the inhomogeneities but can 
nevertheless explicitly calculate the backreaction. 

4 Discussion 

4.1 Potential improvements to the model 

While our model is interesting and internally consistent and perhaps even an improvement 
in sophistication compared to earlier models, what we want is to know is the magnitude and 
nature of the backreaction in the real universe. The most obvious next step in improving our 
model is to add a cosmological constant term, A. With that achieved, it would be possible 
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Figure 7. This figure depicts tiie volume at each gridpoint using the same conventions as figure 6. 
The only difference is that this figure was produced using the model described in section 3.2.4. The 
list of parameters used for this figure also appears in section 3.2.6. See section 3.2.6 for a discussion 
of both figures. 

to do a full Markov Chain Monte Carlo over the parameters of the model to see which set 
of parameters fits the observed data best. If our model, with backreaction, fitted the data 
better than the concordance model, without backreaction, then this would indicate that the 
assumptions we made to solve the various complications were well-motivated. 

However, another complication that is less easily solved would need to be addressed 
before comparisons to the real universe can be meaningfully made. In this paper we have 
calculated the backreaction that arises as structures form over a small range of scales. We 
can vary where this range lies, but we cannot make it arbitrarily wide. To compare the 
model to observation we need a reliable calculation of the total backreaction that arises over 
all scales. In principle this problem can be solved by considering a larger box; however, in 
practical terms, this can never be fully achieved. 

Given the hierarchical nature of structure formation, an ideal solution to this problem 
is a time-dependent grid spacing and box size. The time-dependence could be global, with 
the grid spacing R set as a function of time R{t). R{t) would be chosen based on which 
scales start forming structures at which times and at which times structure growth freezes 
out. Alternatively, the time-dependence could be set locally. There could be a progressive 
"zooming out" as individual gridpoints expand into neighbouring gridpoints. Whenever two 
gridpoints collide it would be possible to stop considering them separately and to recalculate 
the metric from that time onwards as if they were one gridpoint. Precisely what metric to 
use for this combined gridpoint would need to be determined. Either a global or a local time- 
dependence to the grid spacing would face this difficulty. That is, the question of precisely 
what initial conditions, or background model, to use for the zoomed out metric. Should to 
(which was set as a function of h) become a time-dependent quantity? If so, should the h 
used in the transfer function also be dependent on the grid spacing and thus time? 
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If the zooming out problem could be solved it would remove any problems that relate to 
how to fix the over expansion of gridpoints. However, if the zooming out problem is not solved, 
there are some improvements that can be made to our fixing methods. Firstly, it should be 
possible to pick (7i.at fi'om observations, rather than to treat it as a free parameter, as we 
have here. Secondly, we chose to fix our expanding regions to the initial background model. 
This is not ideal. When any region reaches the point where we wish to fix its expansion, 
the average expansion rate, and more importantly, the average volume and density, are not 
the same as those of the initial background, extrapolated to this time. If the universe is 
to reach a new asymptotic FRW state it should be related to the new average density, not 
the one extrapolated from the initial background model. However, doing this is not easy. 
The deviation of the average density from the initial background model changes with time. 
Therefore, the correct background model to fix the expanding regions to also changes with 
time. When the total backreaction is < 10% these issues will be less important than the many 
other features of the model. However, when comparing the model to precise measurements, 
perhaps 10% of 10% is no longer an ignorably small quantity. 

4.2 Conclusions 

Let us close by summarizing our methodology and the conclusions we can draw from it. 
We have approximated the evolving inhomogeneous synchronous gauge metric of a CDM 
$7 = 1 universe with realistic initial conditions using a gradient expansion, and studied the 
backreaction on the average dynamics. This gradient expansion is an approximation to the 
true metric of such a universe which goes beyond the description of standard perturbation 
theory. Instead of being limited by the smallness of the inhomogeneities, it has a temporal 
range of validity set by the initial local 3-curvature. For timescales t > icon the gradient 
expansion does not converge anymore. This does not pose a real problem for studying 
regions with positive initial curvature which collapse at around that timescale and eventually 
virialize. However, it does present limitations in following the late time evolution of under- 
dense regions with negative initial curvature. Nevertheless, if the late time evolution of such 
regions is fixed in some prescribed manner, this approach offers a way to construct a realistic 
model even for the late-time metric of such a universe. 

Our methodology improves over previous approaches in that it is fully relativistic, goes 
beyond perturbation theory and does not treat the universe as a collection of unconnected 
under-dense and over-dense structures with high symmetry. Of course it does require mod- 
eling of the late time behaviour of the under-dense regions which is not naturally captured 
by the gradient expansion but such modeling can be simply parameterized and possibly in- 
formed by observation. We have attempted to achieve this in a simple manner in this paper 
but it can still be made more realistic in a number of ways. 

Given the approximations mentioned, we computed the backreaction of inhomogeneities 
in this model metric. Our qualitative results indicate that backreaction may constitute 
more than a percent effect and is thus highly relevant in future considerations of precision 
cosmology. Of course more realism is needed before definite statements can be made regarding 
the backreaction in our universe. 

There is another aspect in which such a model would aid in understanding the effects 
of inhomogeneities for cosmology. We have focused on the average dynamics but what we 
actually see is light propagating through the inhomogeneous universe and not the average 
scale factor. One could question whether the latter is at all relevant for our observations. 
To really assess the impact of inhomogeneities on observations one should trace light rays 
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through the inhomogeneous spacetime and determine what observers would see [37, 38]. This 
is easy to do in a box described by our metric. In particular we can compute what happens 
to the redshift of photons or the travel time though over-dense and under-dense regions. The 
great advantage of our approach is that this is easily calculable within our model. Therefore, 
even irrespective of whether our model describes the universe entirely accurately, we can 
answer the question of whether a large backreaction in the synchronous gauge corresponds 
to a large shift in the time it takes a photon to traverse the universe and/or the redshift it 
experiences. We will return to this issue and further improvements to our model in future 
work. 
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A The squared metric 

The metric (2.33) is the result of a gradient expansion approximation to the Einstein equa- 
tions for a dust dominated universe. It has however the undesirable feature that this ex- 
pression eventually becomes negative in regions of positive curvature. To overcome this, the 
authors of [28] advocated the following fix: they took the square root of the original metric 
to the desired gradient order and then squared it again. More concretely, in our case this 
procedure yields for the 1st and 2nd order metric 




(A.l) 



7ij 



(A.2) 



where Bij is now 
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Figure 8. The "squared metric" for tlie same configuration as in figure 1. Tlie solid line is the exact 
solution, the dotted line is the zero gradient approximation for the metric (flat FRW or "separate 
universe approximation"), the dashed line is the 1st order metric (2 spatial derivatives) squared and 
the dash-dotted line is the 2nd order metric (4 spatial derivatives) squared. The difference in the 
behavior compared to the straightforward gradient expansion for large times is obvious. 

Note that now there are higher order terms which ensure that the metric is always positive. 

It was argued in [28] that such a procedure improves convergence and in fact can re- 
produce exactly the Szekeres solution studied there. When applied to our spherical test case 
it actually does worse than the straightforward gradient expansion as can be seen from fig- 
ure A. However, its use for studying backreaction has merits, as argued in the main text. 
Furthermore, we see that the (less important for backreaction) collapsing regions are still 
modeled relatively well while the behaviour of the expanding regions at two gradients gets 
regulated since they eventually slow down. In this configuration our recipe for fixing the 
squared metric would result in the collapsing region shrinking down to zero volume (and 
hence be irrelevant for backreaction) while the expanding region would be made to follow the 
background evolution at t ~ 14. This corresponds to a volume ratio between the expanding 
region and the background Kxp/^bg ~ 3.5. 

B Fixing the metric continuously 

To avoid sharp features appearing in Q{ao) (and to a lesser extent Q) it is advantageous 
to smoothly impose the fix of the expanding regions to the background. For the model 
that consists of spherically symmetric regions, introduced in section 3.2.4, the asymptotic 
state that is required is written in eq.(3.13). We choose to to impose the following smooth 
transition 
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imposed for all t > t*, which asymptotically approaches eq.(3.13). Note, 7 = 7(t*,x*). 
Remember that for the model of spherically symmetric regions, 711 = 722 = 733 and 7jj = 
for all i ^ j. Note that at t = t*, ju = 7^/^ and ju = jii{t*) as is required for the transition 
to be smooth. Although (and therefore ci'd) will not be continuous during this transition, 
it is well-defined. By increasing a the transition can be made sharper. This will result in 
sharper features in Q{ao)- By making a smaller the transition can be made to be more 
gradual. However this has the effect of increasing the final constant term in the expression 
for 7jj in eq.(B.l) making the grid point take much longer to reach the asymptotic state. 

For the model introduced in section 3.2.5, the issue of smoothness in the transition is 
not as important. That model uses the squared metric and all of the components in ^^ij. For 
that model, we actually impose the fix to the background at precisely the moment when 7 
at some grid point is the same as the background. That is, the moment when the gridpoint 
stops expanding faster than the background. We still apply the fix described in (B.l), with 
the additional conditions that 7ij = 7ij = for all i / j (and all i > i*). This instantly 
forces the gridpoint into a spherically symmetric state, which is obviously not a smooth or 
even continuous transition; however 7 and 7 do remain continuous and therefore so do ao 
and (Id- 

For the two lower panels of figure 2 we used a = 1. For figure 4 we usd a = 50. For 
figure 5 we used a = 20. For figure 6 we used a = 5. And for figure 7 we used a = \. The 
much larger values of a used for the model of section 3.2.5 is a result of the transition already 
being smooth. That is, using eq.(B.l) instead of eq.(3.13) does not improve the situation. 
Setting a 3> 1 effectively makes the transition instantaneous. However setting a ^ 1 for the 
model of section 3.2.4 makes Q{a£)) much more noisy. 
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